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Chapter  1 
Introduction 


For  as  long  as  quasi-phased  matched  devices  have  been  produced,  there  has  been  a  need  to 
evaluate  and  improve  their  performance.  This  evaluation  has  most  often  been  done  qualitatively 
by  spot  checking  domain  boundaries  of  etched  surfaces  using  a  light  microscope  or 
quantitatively  by  using  the  crystal  in  a  power  measurement  to  determine  conversion  efficiency. 
Spot  checking,  however,  only  gives  coarse  information  on  whether  the  device  will  be  effective, 
and  the  conversion  efficiency  measurement  requires  that  an  experimental  setup  be  available  for 
the  purpose  of  testing  the  crystal.  This  is  not  always  practical,  especially  for  crystal  vendors  who 
serve  a  wide  variety  of  applications  which  would  require  a  host  of  lasers  to  test  the  crystals. 
Furthermore,  as  quasi-phase  matched  designs  become  more  complex,  such  as  using  chirped 
structures^  and  even  aperiodic  structures^,  an  evaluation  of  how  close  the  poled  structure 
reproduces  the  mask  design  is  required. 

A  new  approach  to  better  quantify  the  underlying  domain  structure  was  developed  to  determine 
the  device  performance.  Previously,  the  best  way  to  quantitatively  compare  two  periodically 
poled  crystals  without  testing  with  a  laser  was  by  measuring  numerical  statistics  on  the  domain 
period.^  Typically,  this  would  be  from  a  representative  region  of  the  crystal,  instead  of  using  the 
complete  domain  structure.  In  this  thesis,  the  use  of  image  processing  across  an  entire  z+  or  z- 
surface  of  a  poled  crystal  is  explored  to  characterize  the  domain  structure  and  then  is  compared 
to  the  results  of  a  representative  region.  The  poled  regions  are  mapped  by  finding  regions  of 
maximum  curvature  in  the  image  surface,  or  alternatively  maximum  derivatives,  corresponding 
to  the  domain  boundary  edges  produced  during  the  etching  process,  and  the  mapping  is  used  to 
determine  an  effective  nonlinear  coefficient.  Experimentally  validated  measurements  from  an 
optical  parametric  generation  (OPG)  setup  are  consistent  with  the  calculated  nonlinear 
coefficient  across  the  crystal.  This  technique  would  be  useful  for  evaluating  the  quality  of  device 
iterations  for  either  research  or  production. 

1. 1  Introduction  to  Nonlinear  Optics 

Nonlinear  optical  crystals  have  many  useful  purposes.  One  of  the  most  common  applications  is 
converting  the  wavelength  of  a  laser  to  one  or  more  other  wavelengths.  Indirectly  generating 
wavelengths  of  light  at  atmospheric  transmission  windows,  for  example,  is  useful  for  military 
applications,  such  as  laser  radar,  infrared  countermeasures,  remote  sensing,  and  target 
recognition  and  designation.  Commercial  applications  include  spectroscopy,  displays,  and 
various  medical  technologies. 

Nonlinear  optical  interactions  in  these  crystals  work  by  driving  a  nonlinear  polarization  field,  P, 
often  but  not  necessarily  expressed  as  a  power  series  in  the  field  strength,  E,  shown  in  Eq.  (1)  as 

P(t)  =  +  x'"'E\t)  +  x'''E\t)  +  ...  (1) 

where  %  is  the  susceptibility,  the  tilde  denotes  a  quantity  that  varies  rapidly  in  time,  and  the 
superscript  (n)  denotes  the  n-th  order.  The  first  order  susceptibility  is  also  known  as  the  linear 
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susceptibility  and  is  related  to  the  index  of  refraetion  and  dispersion.  The  foeus  of  this  researeh 
deals  with  processes. 

The  general  idea  of  nonlinear  optical  processes  in  materials  is  that  the  incident  field  strength,  E, 
induces  oscillating  dipole  moments  in  eaeh  of  the  atoms.  A  thorough  deseription  of  the  origin 
and  physics  of  nonlinear  optical  processes  is  given  by  Boyd."^  Due  to  the  nonlinear  nature  of  the 
response,  a  nonlinear  polarization  at  new  frequeneies  is  generated  whieh  ean  radiate  at 
frequeneies  not  present  in  the  incident  radiation  field.  This  coupling  allows  energy  to  be 
transferred  between  different  wavelengths  and  forms  the  basis  of  the  physieal  mechanism  behind 
these  proeesses.  An  isolated  atom  would  radiate  in  the  typical  dipole  radiation  pattern,  but  in  a 
material,  a  large  number  of  dipoles  are  oscillating  with  a  phase  related  to  the  ineident  fields 
producing  a  phased  array  of  dipoles.  With  nonlinear  interaetions,  the  incident  fields  set  up  a 
nonlinear  polarization  which  needs  to  be  phased  up  with  the  freely  propagating  field.  If  the 
dipoles  radiate  in  phase  with  each  other,  then  the  resulting  radiation  will  add  constructively  and 
form  a  narrow  beam.  This  proeess  is  somewhat  analogous  to  a  phased  array  of  antenna  elements 
that  radiate  at  radio  frequencies. 

1.2  Nonlinear  Optical  Interactions 

The  following  derivation  for  the  propagation  of  plane  wave  light  through  a  (lossless)  nonlinear 
optical  material  is  provided  with  extensive  details  by  Boyd."^  Beginning  with  Maxwell’s 
equations,  shown  in  Eq.  (2)  through  Eq.  (5), 

V -0  =  471^  (2) 

V-.8=0  (3) 


VxE 


1  dB 
c  dt 


(4) 


VxH 


1  dD 

c  dt 


(5) 


it  is  assumed  there  are  no  free  eharges  (p  =  0),  no  free  eurrents  (J  =  0),  and  the  material  is 
nonmagnetic  (B  =  H).  The  nonlinear  nature  of  the  material  is  reflected  in  the  displacement  field 
relationship  shown  in  Eq.  (6), 

D  =  E  +  47rP.  (6) 


By  taking  the  curl  of  Eq.  (4)  and  substituting  in  Eq.  (5),  the  equation  eventually  simplifies  for  the 
nonlinear  proeesses  discussed  shortly  into  an  inhomogeneous  wave  equation,  shown  in  Eq.  (7), 


V^E 


d^E 
c"  de 


An  d^P^^ 
dt^ 


where  n  is  the  index  of  refraetion  and  is  the  2"‘^  and  higher  order  polarizations. 


Consider  an  incident  field  made  up  of  two  independent  frequeney  components  producing  a 
second  order  polarization  through  the  nonlinear  susceptibility.  The  incident  field  is  shown  in  Eq. 
(8), 
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E{t)  =  +  E^e-“^‘  +  c.c. ,  (8) 


which  produces  the  polarization,  shown  in  Eq.  (9), 

P(0  =  +c.c.]  (9) 

+2z'^>[e,e; +e,e;1 

These  terms  represent  the  possible  order  nonlinear  interaetions,  but  typieally  only  one  term  in 
Eq.  (9)  will  be  used  sinee  it  is  unlikely  that  more  than  one  term  will  satisfy  the  phase  matching 
condition  simultaneously. 

1.2.1  Difference  Frequency  Generation  (DFG) 

In  one  seeond  order  nonlinear  interaetion  ealled  DEG,  shown  in  Eigure  1,  two  optical  waves  at 
frequeneies  cOp  and  cOj,  representing  the  pump  and  signal  waves,  interaet  in  the  nonlinear  optieal 
material  to  produee  the  idler  wave  at  frequeney  coi  =  cOp  -  coj.  hollowing  the  derivation  in  Boyd, 
ultimately  leads  to  a  set  of  coupled  amplitude  equations,  shown  in  Eq.  (10)  and  Eq.  (11), 


dz 


Smo), 

kc-  ^ 


(10) 


dA. 

dz 


^TncOj  d  ... 


(11) 


where 


E^  =  A 


Ak^z 


m  =  p,  s,  or  i 


Ak  =  kp -k^  -k, 


The  eoupled  amplitude  equation  for  the  pump  wave  is  not  provided  because  the  pump  amplitude 
is  assumed  to  be  undepleted,  or  the  relative  spatial  ehange  is  so  small  as  to  be  negligible. 


Figure  1.  Difference  frequency  generation  a)  system  and  b)  photon  energy  ievel 
diagram. 
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The  phase-matched  condition  is  AA:  =  0  and  any  deviation  reduces  the  efficiency.  The  reason  is 
that  the  exponential  in  Eq.  (11)  containing  Ak  is  an  oscillating  term  in  the  z-direction.  For  DEG 
with  two  undepleted  pumps,  we  can  simply  integrate  this  equation.  The  result  is  that  when  Ak  is 
close  to  zero,  the  spatial  extent  of  the  oscillation  is  larger  than  the  crystal  length  leading  to  a  net 
positive  after  integration.  The  largest  integration  will  occur  when  Ak  =  0,  since  the  exponential 
term  is  equal  to  unity.  Once  the  oscillation’s  period  is  shorter  than  the  crystal  length,  the 
oscillation  tends  to  have  nearly  parts  positive  and  negative,  which  counteract  each  other  so  the 
field  does  not  get  a  chance  to  build  up.  The  result  is  that  for  Ak  away  from  zero,  the  integration 
equals  nearly  zero. 

1.2.2  Optical  Parametric  Generation  (OPG) 

Optical  parametric  generation  is  a  process  where  the  signal  and  idler  fields  are  spontaneously 
generated.  It  is  sometimes  also  known  as  spontaneous  parametric  fluorescence.  The  pump  photon 
still  splits  into  a  signal  and  an  idler  photon,  but  the  signal  and  idler’s  wavelengths  are  determined 
by  the  phase  matching  condition  and  the  energy  conservation  principle.  The  phenomenon  can  be 
described  semi-classically  as  being  started  by  quantum  vacuum  fluctuations  that  allow  emissions 
to  occur  at  any  wavelength  and  in  any  direction  which  satisfy  energy  and  momentum 
conservation.  To  account  for  this,  the  signal  and  idler  can  be  initially  set  to  14  of  a  photon  or  the 
signal  can  be  set  to  one  photon  and  the  idler  set  to  none. 

The  OPG  process  does  not  have  physical  threshold  like  other  optical  processes  since  there  will 
always  be  some  output.  Instead  a  “threshold”  is  determined  by  the  lowest  detectable  output.  The 
pump  threshold,  therefore,  will  be  the  pump  power  (or  energy)  that  produces  a  detectable  output 
(signal  and  idler).  To  calculate  this  threshold,  the  expected  signal  power  collected  by  a  detector 
subtending  an  angle  20  for  a  given  pump  power  assuming  small  gain  is  calculated  from  Eq. 

(12)', 


b\ 


0^ 


(12) 


where 


iTi^Si^n^ripC^ 


and  b  = 


(  dk^ 

f  dk.  ^ 

[d(Dj 

S 

II 

^dco^j 

In  the  equation,  n  is  index  of  refraction,  co  is  angular  frequency,  k  is  wavevector,  c  is  speed  of 
light,  8o  is  the  dielectric  constant  of  free  space,  h  is  Dirac’s  constant,  E  is  length,  P  is  the  power, 
the  subscripts  p,  s,  and  i  represent  the  pump,  signal,  and  idler,  respectively,  and  the  subscript  0 
refers  to  the  frequency  at  perfect  phase  matching.  By  looking  at  the  output  of  the  OPG  at  the 
threshold,  we  ensure  that  the  undepleted  pump  (small  gain)  approximation  is  valid.  When 
translating  laterally  across  a  crystal  using  quasi-phase  matching  (defined  below),  all  of  the 
parameters  will  be  constant  except  for  deff  (due  to  variations  in  the  quality  of  the  quasi-phase 
matching),  so  a  scaled  relative  threshold  is  inversely  proportional  to  dgf/. 
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1.3  Quasi-Phase  Matched  Analysis 

There  are  two  main  ways  to  aehieve  phase  matching:  birefringent  phase  matching  (BPM)  and 
quasi-phase  matching  (QPM).  BPM  uses  the  birefringence  of  a  crystal  to  compensate  for  the 
dispersion  of  the  linear  refractive  indices.  This  typically  involves  orienting  a  crystal  and  laser 
beam  to  achieve  the  desired  interaction.  When  BPM  is  not  normally  possible  or  more  flexibility 
is  desired  (different  output  wavelengths),  QPM  can  also  create  a  desired  phase  matching 
condition.  A  review  of  recent  advances  in  QPM  is  given  by  Hum  and  Fejer.^  The  essential  idea  is 
modulation  of  the  nonlinear  susceptibility  as  a  function  of  propagation  distance  in  a  crystal.  As 
described  with  birefringent  phase  matching  above,  if  Ak  is  not  equal  to  zero,  the  field  will  build 
up  and  then  decay.  The  reason  for  this  is  that  the  nonlinear  polarization  at  a  given  position  in  the 
crystal  radiates  and  this  radiation  can  interfere  with  the  freely  propagating  field  generated  earlier 
in  the  crystal.  Hence  the  field  builds  up  for  one  coherence  length  and  then  decays  for  another 
coherence  length.  At  the  point  where  the  field  starts  to  decay,  one  can  flip  the  nonlinear 
polarization  to  essentially  reset  the  phasing  so  that  the  field  continues  to  build  up.  For  a  sign 
inversion  of  the  d  coefficient,  this  results  in  a  180  degree  phase  shift  of  the  nonlinear 
polarization.  By  spacing  these  inversions  by  one  coherence  length,  the  field  will  grow,  although 
at  a  slower  rate  than  perfect  phase  matching.  The  size,  shape,  and  orientation  of  the  designed 
crystal  structure  determine  which  output  wavelengths  satisfy  the  phase  matching  condition, 
unlike  BPM,  where  optical  material  properties  limit  the  possibilities.  Also,  larger  nonlinear 
coefficients,  unavailable  for  BPM,  can  be  selected. 

QPM  materials  can  be  produced  in  several  ways.  Originally  devised  by  Armstrong  et  al.^,  early 
attempts  at  producing  QPM  materials  involved  stacking  rotated  crystals  and  diffusion  bonding 
them  together.  This  idea  was  abandoned  for  much  better  results  using  lithography  to 
periodically  pole  ferroelectrics  with  a  large  electric  field.  The  periodical  poling  of  a  crystal  is 
the  most  common  method  for  creating  a  QPM  material.  In  periodic  poling,  an  original  bulk 
crystal  with  the  z+  axis  facing  upward  is  transformed  into  a  crystal  with  domains,  or  region  with 
a  uniform  crystal  orientation,  that  alternate  between  z+  axis  being  up  and  down,  shown  in  Figure 
2.  The  effect  of  this  transformation  is  to  produce  domains  that  alternate  between  +deff  and  -deff 
across  the  crystal,  effectively  converting  the  d  coefficient  from  a  constant  to  a  spatially  varying 
pattern. 


Figure  2.  Representation  of  QPM  through  electric  field  periodic  poling.  A) 
Homogeneous  single  crystal  (before  poling)  and  b)  a  periodically  poled  material 
with  alternating  inverted  crystal  orientations  and  A  spatial  period. 
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The  effect  of  a  quasi-phase  matched  material  on  the  DFG  and  OPG  processes  is  to  change  the 
phase  matching  condition  and  the  relative  strength  of  the  d  coefficient.  The  first  step  is  to  include 
a  spatially  varying  d  coefficient,  d(z)  =  g(z)deff,  in  the  coupled  amplitude  equations,  Eq.  (10)  and 
Eq.  (11),  to  form  Eq.  (13)  and  Eq.  (14), 

^  =  (13) 

dz 

^  =  7,d„A:g{z)e-“--  (14) 


where 


^  ^TUCO^ 

AA:'=  — AA:  =  A:  +  k,  —  k. 


m  =  s,  1 


Eor  DEG  with  a  strong  undepleted  pump  and  signal  (slowly  varying  amplitude  approximation), 
the  idler  amplitude  as  it  leaves  the  crystal  simplifies  to  several  constants  multiplied  by  the 
Eourier  transform  of  the  g(z),  as  described  by  Eejer  et  al.  and  as  shown  in  Eq.  (15), 


4(i)  =  r,<,,j;jg(z)e-“'--& 


=  3{g(z)} 


(15) 


The  limits  on  the  integral  result  from  the  finite  extent  of  the  nonlinear  coefficient  (within  the 
crystal).  The  result  from  normal  phase  matching  when  g(z)  =  1  within  the  crystal  and  g(z)  =  0 
everywhere  else  is  shown  in  Eq.  (16), 

4(L)  =  A^'E/2).  (16) 

To  solve  for  a  spatially  periodic  nonlinear  coefficient,  g(z)  is  first  written  in  terms  of  its  Eourier 
series,  shown  in  Eq.  (17), 

g(z)  =  XG.e“-'-  (17) 


where  K^=ln  !  K  is  the  fundamental  spatial  frequency,  and  A  is  the  spatial  period.  The  ideal 
structure  for  periodic  poling  would  be  a  square  wave  modulation  from  -1  to  1  with  a  50%  duty 
cycle.  In  this  case,  the  Gn=  2/(n7T:)  and  the  new  coupled  amplitude  equation  looks  like  Eq.  (18), 

,  ”  ■  .  (18) 

0  « 
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where  AAT  =  k^+  -  k^  -K^,  indieating  a  new  phase  matehing  eondition.  Only  one  of  the  Gn 

eoeffieients  will  be  close  enough  to  AK  =  0  to  dominate  the  integral,  so  the  series  collapses  to  a 
single  value  based  on  the  desired  interaction.  The  final  result  is  similar  to  Eq.  (16),  as  seen  in  Eq. 
(19), 


A,{L)  =  r,d„GXl 

0 


e-‘^^dz 


=  smc(AKL/2) 


(19) 


where  dn  =  deffGn. 

A  non-ideal  periodic  structure  such  as  square  wave  with  duty  cycle  variations  would  have  a 
similar  Eourier  spectrum  to  the  ideal  square  wave,  except  each  ideal  peak  (corresponding  to  each 
Gn)  would  have  a  reduced  magnitude  and  a  broadened  line.  These  can  be  calculated  by  changing 
the  spatial  period  to  the  crystal  length  for  the  Eourier  series  calculation.  If  the  variation  in  duty 
cycle  is  modest,  the  broadening  would  be  small  compared  to  the  dominant  peak  and  the  phase 
matching  bandwidth  (sine  function)  would  still  be  small.  Any  nonlinear  coefficient  contribution 
from  off  the  AK=0  peak  would  diminish  quickly,  allowing  the  main  peak  to  remain  dominant. 
Even  with  a  non-ideal  square  wave,  the  amplitude  of  the  main  peak  of  a  desired  interaction  will 
be  used  for  measuring  its  nonlinear  coefficient  strength.  Eor  the  same  reasons,  g(z)  in  the  OPG 
coupled  amplitude  equations  can  be  replaced  with  the  value  of  the  main  peak  near  phase 
matching.  Eor  more  complex  patterns  like  chirped  structures,  g(z)  would  no  longer  be  able  to  be 
reduced  to  a  single  number  and  a  more  sophisticated  calculation  would  be  required  to  determine 
the  relative  efficiency  and  bandwidth. 

1.4  Visualization  Techniques 

The  quantitative  characterization  of  QPM  nonlinear  optical  materials  requires  the  use  of  a 
visualization  technique  before  it  can  be  digitized  and  processed  to  calculate  performance  metrics. 
A  review  of  many  of  the  visualization  techniques  for  ferroelectric  domains  is  given  by  Soergel. 
Many  of  those  methods  of  visualization  were  considered  but  either  did  not  have  enough 
resolution  (better  than  1  qm)  or  the  acquisition  of  an  entire  crystal  would  have  taken  too  long. 

The  options  examined  were  using  a  SEM,  a  new  technique  involving  doping  Rb  into  KTP,  and 
visual  light  microscopy  of  differentially  etched  Eithium  niobate  and  KTP  crystals. 

1.4.1  Scanning  Electron  Microscope  (SEM) 

SEM  imaging  is  one  domain  imaging  technique  that  provides  excellent  resolution  and  is 
relatively  easy  and  fast  to  perform  for  crystals  such  as  lithium  tantalate  and  KTP  .  Elsing  a 
SEM  with  certain  settings  produces  a  domain  contrast  due  to  screening  of  charge  which  differs 
between  the  z+  surface  and  the  z-  surface  due  to  heating  from  the  e-beam  and  the  pyroelectric 
effect  as  seen  in  Eigure  3.  With  other  settings,  the  domain  walls  become  dark  or  bright  due  to  a 
change  in  surface  profile  height  from  the  piezoelectric  effect  and  the  orientation  of  the  detector 
as  seen  in  Eigure  4.  It  is  possible  to  image  an  entire  crystal  using  either  method. 
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Figure  3.  SEM  image  of  periodicaliy  poied  KTP  showing  intensity  contrast. 


Figure  4.  Close-up  SEM  image  of  periodically  poled  KTP  showing  boundary  contrast. 

1.4.2  Rb-dopedKTP 

A  new  technique  for  domain  visualization  was  developed  for  KTP.  The  idea  was  to  dope  KTP 
with  rubidium  by  dipping  it  into  a  RbNOs  melt,  which  differentially  dopes  the  surface  (up  to  a 
micron  in  depth)  depending  on  the  domain  (z+  or  z-)  which  is  exposed.  By  taking  a  poled  KTP 
sample,  doping  it,  and  then  examining  the  surface  for  its  level  of  Rb  concentration,  a  map  of  the 
domain  structure  can  be  produced.  The  surface  composition  can  be  measured  using  a  SEM 
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incorporating  an  energy  dispersive  spectroscopy  (EDS)  system.  This  system  performs  x-ray 
microanalysis  which  can  detect  the  various  elements  within  a  material.  The  results  were  not  as 
impressive  as  hoped,  however.  Several  immersion  times  were  tested  with  different  pieces  of  a 
commercially  bought  periodically  poled  KTP  sample.  A  100  nm  layer  of  gold  was  sputtered  on 
the  crystal  pieces  before  viewing  under  the  SEM  to  avoid  the  accumulation  of  static  electric 
charge  during  electron  irradiation.  The  thickness  of  the  gold  was  smaller  than  the  depth  of  the 
electron  beam’s  interaction  volume  within  the  sample,  which  can  extend  up  to  5  pm,  so  all  of  the 
elements  within  the  compound  were  still  detectable.  The  resulting  images  in  all  cases  were 
similar,  with  the  best  regular  SEM  vs.  Rb  image  comparison  shown  in  Eigure  5.  Not  only  did  the 
Rb  images  have  poor  contrast  and  definition,  but  the  image  acquisition  took  much  longer  than 
anticipated  (hours)  for  an  image  of  a  couple  hundred  microns  on  a  side.  This  was  due  to  the  low 
detector  count  rate  for  the  energy  level  corresponding  to  Rb. 
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Figure  5.  Side-by-side  comparison  of  reguiar  SEM  and  EDS  image  of  Rb  density. 
1.4.3  Chemical  Etch 

The  differential  etching  of  periodically  poled  lithium  niobate  and  KTP  allows  the  domains  of 
these  crystals  to  be  seen  easily  using  a  visible  light  microscope.  The  acid  HP  etches  lithium 
niobate  crystal  surfaces  at  different  rates  for  the  z+  and  z-  surfaces.  The  discontinuities  in  the 
surface  profile  produce  dark  edges  under  a  microscope  as  seen  in  Eigure  6.  A  mixture  of  KNO3 
and  KOH  in  heated  deionized  water  etches  the  surfaces  of  KTP  at  different  rate  for  the  z+  and  z- 
surfaces.  A  discontinuity  is  visible  for  etched  KTP,  also,  as  seen  in  Eigure  7.  Although  these 
methods  damage  the  surfaces  of  the  crystals  to  some  degree,  the  etching  process  is  quick,  the 
resolution  is  good,  and  taking  pictures  of  an  entire  crystal  under  microscope  is  feasible.  This 
technique  was  chosen  for  imaging  the  QPM  domains. 
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Figure  6.  Picture  of  etched  periodically  poled  lithium  niobate  under  a  visibie  iight 
microscope. 


Figure  7.  Picture  of  etched  periodicaiiy  poied  KTP  under  a  visible  light 
microscope. 
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Chapter  2 
Methodology 

The  main  tasks  for  this  experiment  are  (a)  to  obtain  an  image  of  the  entire  crystal,  (b)  to  map  the 
crystal  domains  by  determining  the  positions  of  the  domain  boundaries  through  image 
processing,  (c)  to  calculate  the  modified  effective  nonlinear  coefficient,  dgff,  and  then  (d)  to 
compare  it  to  measurements  based  on  an  OPG  setup.  The  relative  spatial  variations  in  threshold 
and  output  (signal  and  idler)  energy  measurements  should  match  the  calculations  based  on  the 
calculated  deff. 

2.1  Image  Registration 

To  map  the  crystal  domains  of  a  periodically  poled  material,  pictures  of  an  entire  z+  or  z-  face 
must  be  acquired.  This  can  be  accomplished  by  taking  a  series  of  slightly  overlapping  pictures 
(using  a  motorized  stage  or  manually)  of  an  etched  material,  in  this  case  periodically  poled 
lithium  niobate  (PPLN),  using  a  reasonably  high  magnification  optical  microscope  (Nikon  AZ- 
100).  Some  software  packages  will  combine  this  into  a  single  image  resulting  in  a  single  massive 
file,  but  such  a  large  file  can  be  difficult  to  use.  An  alternative  is  to  calculate  how  to  align  the 
pictures  by  registering  them  together.  There  are  many  strategies  for  registration  ,  but  the  correct 
technique  often  depends  on  the  circumstances.  Since  poled  materials  are  often  periodic,  the  use 
of  cross-correlation  alone  can  be  complicated  due  to  the  many  peaks  that  would  result. 

A  feature-based  registration  technique  is  more  appropriate  for  our  application  since  poling 
imperfections  and  dirt  in  the  overlapping  regions  can  be  used  as  features  for  aligning.  The 
features  can  be  extracted  manually  or  by  looking  for  locations  where  there  is  a  maximum 
curvature  or  derivative  along  both  coordinate  axes  to  find  a  sharp  comer  in  the  domain  boundary 
or  the  location  of  a  speck  of  dust.  Since  the  images  are  obtained  using  a  raster  pattern  with  a 
consistent  amount  of  overlap,  the  translation  relationship  between  the  features  can  be  estimated. 
The  estimate  will  be  more  accurate  if  the  images  are  collected  using  a  motorized  stage  with  a 
specified  step  size  due  to  its  accuracy  and  precision.  In  this  initial  investigation,  manual 
registration  was  used  before  any  image  processing  was  performed. 

Before  the  images  can  be  combined,  any  distortion  or  aberrations  should  be  corrected,  if 
necessary,  to  avoid  an  incorrect  registration.  While  most  aberrations  lead  to  an  image  that 
appears  fuzzy  or  out  of  focus,  our  images  were  clear  and  in  focus.  It  is  possible  to  see  a  small 
amount  of  lateral  chromatic  aberration,  but  this  can  be  avoided  by  only  using  the  green  color 
component  and  calibrating  the  vertical  and  horizontal  scales  accordingly.  Scale  of  the  images 
was  not  important  for  the  measurements  in  this  experiment,  so  it  was  not  determined.  Barrel  or 
pincushion  distortion  would  likely  be  the  only  other  significant  issue,  but  it  was  not  noticeable  in 
the  many  nearly  straight  lines  of  the  domains.  Also,  the  amount  of  this  distortion  was  calculated 
using  a  calibration  slide,  which  has  many  evenly  spaced,  straight  lines.  The  amount  of  distortion 
(percent  difference  between  actual  and  expected  spacing  of  lines)  did  not  exceed  0.5%,  which 
was  believed  to  be  a  negligible  factor  at  this  point  and  therefore  was  not  corrected.  The  distortion 
would  be  more  of  a  concern  when  the  wavelengths  involved  are  important,  which  was  not  the 
case  in  this  experiment.  Distortion  will  be  examined  more  closely  as  the  development  of  image 
processing  techniques  progresses. 
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By  simply  scanning  the  microscope  across  the  crystal,  the  transformation  to  relate  the  positions 
in  one  image  to  another  should  only  involve  translation.  There  would  not  be  any  scaling  and  the 
rotation  should  be  negligible.  However,  this  will  depend  on  the  translation  stage,  and  care  should 
be  taken  to  make  sure  the  stage  has  minimal  piteh  or  yaw  variations.  After  estimating  the 
adjaeent  images  position,  a  eross-eorrelation  using  a  small  window  and  centered  at  eaeh  feature 
was  performed  against  its  corresponding  location  in  the  other  image.  The  resulting  peak  provided 
the  relative  location  from  the  initial  guess.  The  translation  was  determined  for  all  images  relative 
to  the  first  image.  The  registration  eould  be  further  refined  by  using  the  domain  boundaries  found 
later  as  features  to  mateh,  but  this  was  not  done  in  this  present  work. 

2.2  Image  Processing 

Using  a  representative  mieroscope  image  of  the  etehed  lithium  niobate,  shown  in  Figure  8,  the 
image  proeessing  techniques  are  deseribed  in  detail  below.  A  fdtering  teehnique  using  a  sampled 
radial  Gaussian  function,  G,  shown  in  Eq.  (20), 

=  (20) 

Int 

where  x  and  y  are  integers  and  t  is  a  thickness  parameter,  helps  to  smooth  out  noise  without 
adding  any  artifacts  to  the  image.  By  convolving  this  function  with  the  image,  I,  shown  in  Eq. 
(21), 

L{x,y,t)  =  G{x,y,t)*  I{x,y) ,  (21) 

the  result,  known  as  a  scale-space  representation  of  the  image,  L,  is  obtained.  A  scale-space 
representation  is  useful  as  a  preproeessing  step  because  it  smoothes  features  smaller  than  the 
seale  t,  so  that  only  details  of  a  certain  scale  are  seen.  Usually,  multiple  scales  of  the  seale-space 
representation  are  used  to  differentiate  features  by  their  relative  size  in  images.  In  our  case,  since 
the  lines  are  a  consistent  size,  we  will  only  use  one  of  the  seale-spaee  representations,  although 
multiple  seales  eould  be  used  in  the  future  to  exclude  large  features  from  eonsideration  or 
establish  the  crystal’s  outer  edges. 

The  visual  effeet  of  this  type  of  filter  is  seen  in  Eigure  9.  The  effect  of  the  Gaussian  filter  is  also 
seen  in  line  scans  shown  in  Eigure  10.  The  larger  Gaussian  filter  bloeks  are  mueh  slower  due  to 
the  convolution  operation.  An  alternative  to  the  convolution  is  to  Eourier  transform  the  image 
and  multiply  by  the  Eourier  transform  of  the  Gaussian  funetion,  whieh  is  another  Gaussian 
function.  Eor  this  work,  the  7x7  bloek  filter  (t  =  2)  was  used  because  it  reduces  much  of  the  noise 
while  retaining  the  overall  shape  and  amplitude  of  the  peaks. 
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Figure  8.  Etched  lithium  niobate  test  sample  image. 


Figure  9.  Section  of  original  image  (left),  after  t  =  2,  7x7  pixel  Gaussian  filter 
(middle),  and  after  t  =  10,  61x61  pixel  Gaussian  filter  (right). 
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Figure  10.  Comparison  of  iine  scans  (row  406)  of  red  component  after  Gaussian 
fi  iters. 

The  next  step  after  the  Gaussian  filter  is  to  extraet  useful  information  from  the  image  using 
partial  derivatives  to  find  ridges  or  valleys  representing  domain  boundaries.  A  ridge  or  valley  is  a 
local  maximum  or  minimum  in  at  least  one  dimension,  which  can  be  visualized  as  an  elongated 
object  at  a  different  intensity  level  than  pixels  around  it.  Ridge  and  valley  detection  is  commonly 
used  to  map  out  roads  or  water  systems  in  aerial  photographs  or  to  locate  blood  vessels  in 
medical  images.  There  are  several  techniques  that  use  partial  derivatives  that  are  used  to  find 
ridges  and  valleys.  One  way  to  accomplish  this  is  to  simply  use  the  partial  derivatives  along  the 
horizontal  and  vertical  directions  (Lx  and  Ly)  of  the  image.  Another  option  is  to  calculate  partial 
derivatives  aligned  with  the  principal  curvatures  at  each  pixel  of  a  height  map  of  the  image 
(intensity  is  the  height).  The  key  is  to  transform  the  image  coordinate  system  (x,y)  into  the  p- 
and  q-  directions  of  the  principal  curvatures  by  rotating  by  an  angle  P  defined  by  Eq.  (22)  and 
Eq.  (23) 


COSy^  = 


sin  P  = 
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A  characteristic  of  this  coordinate  system  is  that  Lpq  =  0.  The  assoeiated  direetional  derivative 
operators  are  given  by  Eq.  (24)  and  Eq.  (25), 


dp  =sm{/3)d^-cos{/3)dy. 

(24) 

S,  =cos{fi)d^+sm{fi)d^. 

(25) 

To  deteet  a  ridge  or  valley  with  this  coordinate  system,  a  set  of  conditions. 


Ridge 


Ypp  I  -  I 


or 


Valley  \ 


The  (x,y)  and  the  (p,q)  coordinate  systems  produce  nearly  equivalent  results.  While  the  (p,q) 
eoordinate  system  was  used  for  image  proeessing  of  the  whole  crystal,  ridge  and  valley  detection 
will  be  demonstrated  below  for  a  single  image  using  the  (x,y)  coordinate  system.  The  proeedure 
is  essentially  the  same. 

An  additional  step  of  applying  a  5x5  wiener  filter  after  calculating  each  derivative  helps 
minimize  noise.  A  wiener  image  filter  is  a  low-pass  adaptive  filter  that  using  statisties  from  the 
underlying  image  to  estimate  the  noise  in  the  image,  in  order  to  filter  it  out.  An  original  section 
of  the  image  and  the  images  of  the  derivatives  before  and  after  the  wiener  filter  are  shown  in 
Eigure  12  through  Eigure  19. 


Figure  11.  Original  section  of  image. 
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Figure  12.  Lx  before  wiener  fiiter. 


Figure  14.  Ly  before  wiener  fiiter. 


Figure  15.  Ly  after  wiener  fiiter. 


Figure  18.  Lyy  before  wiener  filter. 


Figure  13.  Lx  after  wiener  filter. 


Figure  16.  Lxx  before  wiener  filter. 


Figure  17.  Lxx  after  wiener  filter. 


Figure  19.  Lyy  after  wiener  filter. 


Assuming  the  domain  boundaries  are  mostly  vertieal,  the  easiest  way  to  determine  the  position  of 
the  boundary  is  to  sean  line  by  line  in  the  image,  looking  for  the  zeros  in  Lx.  To  avoid  small 
variations  and  noise,  thresholds  for  the  first  and  seeond  derivative  must  be  exceeded  within  a 
window  around  a  point  of  interest  stretching  several  pixels.  To  ignore  points  that  are  not  on  a 
line,  a  measure  of  probability  of  a  blob  (such  as  dirt)  was  calculated  using  Eq.  (26), 

Blob  =  L^^*L^,  (26) 

where  Lxx  exceeds  a  threshold,  otherwise  the  value  is  0.  An  image  of  the  blob  value  after  wiener 
filtering  is  shown  in  Figure  20  for  the  same  region  as  in  Figure  12  through  Figure  19. 


Figure  20.  Image  of  the  blob  measurement 

The  potential  ridge  and  valley  points  overlaid  on  the  original  image  are  shown  in  Figure  21 . 

To  further  reduce  potential  ridge  and  valley  points  to  only  include  points  on  the  domain 
boundaries,  a  line  detector  was  constructed.  This  process  determines  whether  a  point  appeared  to 
be  part  of  a  line.  First,  a  small  subsection  of  the  image  around  the  point  was  compared  to  the 
statistical  mode  of  that  subsection  multiplied  by  a  factor  that  was  manually  adjusted  for  the  best 
results.  A  threshold  image  was  created  and  regions  of  contiguous  pixels  around  the  border  that 
exceeded  the  threshold  were  counted  as  long  as  they  were  not  too  small  or  too  large  (parameters 
tuned  for  the  best  results).  A  line  would  have  two  regions  at  the  border  where  the  line  crossed 
into  and  out  of  the  subsection.  The  line  did  not  have  to  be  straight  to  satisfy  this  condition.  The 
only  other  condition  was  that  the  number  of  pixels  above  threshold  did  not  exceed  a  certain 
percentage  of  the  total  pixels  in  the  subsection.  Again,  the  percentage  was  manually  tuned  for  the 
best  results.  This  made  it  unlikely  that  the  pixels  in  the  middle  of  the  threshold  image  showed 
anything  other  than  a  line.  The  procedures  did  not  guarantee  that  the  point  was  part  of  a  line,  but 
it  established  a  good  candidate  that  excluded  uncertain  features.  Usually,  only  a  few  points 
would  be  removed  from  the  list  of  potential  points,  but  these  were  away  from  the  domain 
boundaries,  which  can  be  confusing  for  the  next  few  processes. 

The  reduced  set  of  points  was  next  organized  into  lines  by  scanning  several  rows  in  the  middle  of 
the  image  to  determine  the  approximate  number  lines.  This  was  accomplished  by  adding  up  the 
points  within  groups  of  3  columns  and  looking  for  clusters  of  points  above  a  threshold.  Each 
cluster  was  typically  separate  for  each  domain  boundary.  This  would  not  necessary  work  every 
time,  but  in  this  experiment  it  did  because  the  boundary  lines  are  straight,  vertical,  and  relatively 
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thin  compared  to  the  domain  period.  The  algorithms  will  need  to  be  more  robust  as  they  are 
further  developed. 

Next,  a  point  within  each  cluster  was  connected  to  other  points  by  searehing  up  and  down  within 
a  window  of  a  few  pixels  to  add  more  points  to  that  line.  Here,  the  assumption  is  that  the 
boundary  points  should  be  connected  because  they  are  near  to  each  other.  After  the  lines  were 
established,  linear  interpolation  was  used  to  fill  in  the  gaps. 

By  filling  in  regions  between  the  boundary  lines  with  the  correct  color  (black  or  white),  an 
alternating  binary  map  of  the  domain  structure  for  a  single  image  is  created,  as  shown  in  Figure 
22.  Black  was  for  a  negative  d  coeffieient  and  white  was  for  a  positive  coefficient.  Using  row 
1000  as  an  example,  the  Fourier  transform  of  the  binary  map  is  shown  in  Figure  23,  for  different 
levels  of  zero-padding.  Zero-padding  is  adding  extra  zeros  at  the  end  of  the  array  of  numbers  to 
be  transformed,  which  increases  the  resolution  in  spatial  frequency  space.  The  higher  amounts  of 
zero-padding  give  the  most  accurate  value  of  the  lowest  frequency  peak  in  the  Fourier  transform. 
The  peak  value  can  then  be  used  in  the  OPG  threshold  calculation  to  find  the  relative  pump 
threshold. 


Figure  21.  Boundary  lines  with  linear  interpolation  shown  as  red  dashed  line. 
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Figure  22.  Binary  map  of  domain  structure. 


Figure  23.  Fourier  transform  of  map  at  row  1000  using  various  levels  of  zero- 
padding. 

When  connecting  the  images  together,  the  complete  lines  from  each  image  were  combined  one  at 
a  time  into  a  comprehensive  list.  When  each  new  line  was  added,  it  was  compared  to  all  of  the 
lines  currently  in  the  list  to  see  if  it  overlapped  any  by  calculating  a  mean  difference  between 
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points  on  each  pair  of  lines.  The  registration  information  was  added  to  the  points  making  up  the 
lines  first  and  then  they  were  linearly  interpolated  to  integer  pixel  locations  for  direct 
comparison.  If  the  mean  difference  quantity  was  less  than  a  threshold  (a  value  less  than  1/6  of  the 
domain  period  was  chosen),  then  the  pair  of  lines  were  considered  equivalent  and  the  latest  line 
was  not  included  in  the  comprehensive  list.  Each  row  of  images  was  processed  in  this  way  until 
all  the  domain  boundary  lines  were  included.  The  same  process  for  creating  a  binary  map  and 
calculating  the  EFT  of  a  single  image  was  then  extended  to  a  whole  row  of  images. 

There  were  several  thresholds  and  parameters  that  need  to  be  determined  for  these  procedures  to 
work  correctly.  The  best  way  to  accomplish  this  would  be  to  start  with  a  single  image  and 
manually  adjust  the  parameters  until  each  process  that  uses  the  parameters  produces  the  expected 
results.  The  parameters  should  easily  generalize  to  the  whole  batch  of  images,  assuming  they 
were  all  acquired  similarly.  Eventually,  one  could  conceivably  create  an  algorithm  to 
automatically  tune  the  parameters  by  knowing  the  scale  and  the  design  of  the  device,  and  then 
comparing  the  results  at  each  step  until  they  best  match  what  would  be  expected  for  a  given 
design. 

2.3  OPG  Setup 

An  OPG  setup  was  used  to  validate  the  accuracy  of  the  image  processing  routines.  The  setup  is 
shown  in  Figure  24.  Where  the  1.064|um  wavelength  laser  beam  reaches  the  PPLN  crystal,  the 
beam  waist  (1/e^  diameter)  was  measured  to  be  approximately  0.4  mm.  After  aligning  the  crystal 
for  efficient  conversion  by  ensuring  that  the  back  reflections  were  collinear  with  the  incident 
beam,  the  signal  and  idler  were  directed  to  an  energy  meter  for  measurement.  At  this  orientation, 
the  Fresnel  reflections  from  the  crystal  facets  lead  to  a  low-finesse  monolithic  OPO  which 
slightly  lowers  the  threshold  as  compared  to  when  the  crystal  is  rotated  away  from  the  resonant 
condition.  To  determine  the  threshold,  the  pump  energy  was  adjusted  until  the  combined  signal 
and  idler  energy  reached  1  qJ.  The  pump  energy  was  measured  when  that  condition  was  met  and 
then  the  whole  process  was  repeated  several  times  while  translating  the  crystal  laterally. 


Lens  Idler  beams 

Figure  24.  OPG  experiment  setup. 

2.4  Calculated  and  measured  OPG  threshold 

The  crystal’s  performance  was  compared  to  the  prediction  based  on  mapping  domain  boundaries 
of  the  entire  crystal  top  surface  with  the  image  processing  techniques  described  previously.  The 
measured  and  calculated  pump  threshold  across  the  width  of  the  crystal  is  shown  in  Figure  25. 
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Since  the  calculation  was  relative,  the  constant  of  proportionality  was  adjusted  until  the  scale  of 
the  calculated  threshold  matched  the  measured  thresholds  visually. 


Comparison  of  experimental  and  calculated  pump  threshold 


Figure  25.  Calculated  and  measured  threshold  of  PPLN  crystal. 

For  the  current  experiment,  the  mosaic  of  the  characterized  crystal  was  comprised  of  38  columns 
and  3  rows  of  individual  microscope  images.  The  columns  were  registered  together,  but  for  this 
thesis,  the  rows  have  not  yet  been  connected,  leaving  gaps  in  the  calculated  threshold  data  where 
the  domain  map  is  split  between  rows  in  the  mosaic.  The  reason  for  this  is  that  successive  images 
were  overlapped  significantly  in  the  direction  parallel  to  the  gratings,  but  also  had  a  slight 
translation  perpendicular  to  this  direction,  a  situation  demonstrated  in  Figure  26.  Once  the  rows 
are  connected  together,  the  calculated  threshold  data  will  be  continuous  across  the  crystal.  The 
calculations  show  a  result  consistent  with  the  performance  found  with  the  OPG  setup.  This 
important  proof  of  concept  shows  promise  for  performing  full-scale  realistic  simulations  of 
quasi-phase  matched  nonlinear  crystals. 


21 


Column  4 


Column  2 


Column  3 


Figure  26.  Representative  diagram  of  the  mosaic  of  images. 
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Chapter  3 
Conclusions 

A  technique  for  quantitatively  evaluating  the  poling  quality  of  quasi-phase  matched  crystals  has 
been  demonstrated.  After  stitching  a  mosaic  of  microscope  images  together,  the  locations  of 
crystal  domains  were  accurately  determined  through  a  series  of  imaging  processing  steps.  The 
steps  include  registering  the  images,  Gaussian  fdtering,  using  derivates  to  find  ridges  or  valleys, 
isolating  those  potential  ridges  or  valleys  into  lines,  connected  the  lines,  merging  the  lines  from 
all  the  images,  and  then  creating  a  map  from  these  boundary  lines.  The  full  crystal  map  was  used 
to  calculate  a  relative  dgff  and  then  subsequently  was  compared  experimentally  to  an  OPG 
threshold  measurement  as  a  function  of  position  in  the  crystal.  The  mapping  of  dgff  to  threshold 
matched  the  OPG  measurements,  validating  this  method  as  a  useful  means  to  evaluate  nonlinear 
optical  crystal  performance.  Beyond  this  proof  of  concept,  the  data  extracted  from  the  images 
could  easily  be  applied  to  other  processes  like  second  harmonic  generation  and  be  used  for  any 
laser  wavelength. 

There  are  several  avenues  for  future  research.  In  our  case,  the  crystal  was  selectively  etched 
periodically  poled  lithium  niobate,  but  it  may  be  possible  to  extend  this  technique  to  other 
materials.  For  example,  there  should  be  little  problem  applying  this  to  lithium  tantalate  due  to  its 
similarity.  Also,  a  properly  prepared  orientation-patterned  GaAs  (side  polished)  would  be  a 
possibility.  Alternative  modes  of  obtaining  the  images  are  also  feasible.  Under  certain  viewing 
conditions,  an  SEM  is  capable  of  producing  usable  images  for  mapping  domain  boundaries. 

Once  the  structure  has  been  mapped,  realistic  nonlinear  crystals  can  be  modeled,  supporting 
complex  beam  propagation  simulations. 
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Appendix:  Matlab  Code 


Listing  1:  imProcess.m 

%  Runs  image  processing  routines. 

g, 

0 

%  After  collecting  the  images,  each  image  is  first  Gaussian 
%  filtered,  and  then  all  the  partial  derivatives  and  the 
%  blob  function  are  calculated.  These  are  used  to  find  the 
%  potential  boundary  points,  some  of  which  are  eliminated  because 
%  they  don't  pass  the  IsLine  test.  The  points  are  collected  into 
%  lines  and  then  are  displayed  on  the  images.  All  the  data  for  each 
%  image  is  saved. 


%  Initialize  parameters 

umperpx  =  0.14244. *2; 

w  =  2 ; 

th  =  1  ; 

th2  =  1; 

thb  =  0.25; 

rows  =  3 ; 


Scale  factor  microns/pixel 
Window  for  FindBoundary 
Threshold  for  1st  deriv. 
Threshold  for  2nd  deriv. 
Threshold  for  blob 
Rows  in  superimage 
Columns  in  superimage 


cols  =  38; 

%  Gets  all  the  filenames  for  the  images 
a  =  dir ( ' * . bmp ' ) ; 
for  cnt  =  l:size(a,l) 

fn (floor ( (cnt-1) /cols) +1, mod (cnt-1, cols) +1, : )  =  a (cnt, 1 ) . name; 

end 


%  Iterate  through  all  columns  in  a  row 
irow  =  1; 
for  icol  =  l:cols 
%  Clear  old  variables 

clear  blob  iml  imtestl  Lp  Lpp  Lq  Lqq  linedown  linep  lines  lines2  linesz 
iineup  posx  posx2  pxsz  px2sz  tposx  tposx2  tpx2sz  tpxsz  x  y 

tic 

%  Reading  in  image  and  applying  initial  filters 

disp ([' Reading  '  reshape ( fn ( irow, icol ,:), 1 , 1 6 )  '  and  applying  filters']) 

iml  =  imread ( fn ( irow, icol ,:)) ;  iml  =  iml(:,:,2); 

[imtestl  Lp  Lq  Lpp  Lqq  blob]  =  Filter(iml); 

%  Find  the  potential  boundary  points 

[posx  posx2  pxsz  px2sz]  =  FindBoundary (Lp, Lq, Lpp, Lqq, blob, w, th, th2 , thb) ; 

disp ( ' Done ' ) ; 

toe 


%  Gets  rid  of  points  that  do  not  pass  the  IsLine  test 
disp ( ' Eiiminating  points  that  are  not  lines'); 

%  Ridges  first 
tposx  =  [ ] ; 

tpxsz  =  zeros (1, size (imtestl, 1) ) ; 
for  i  =  1 : size (imtestl, 1) 
for  j  =  1 rpxsz (i) 

if  (IsLine (imtestl (max (1,  (i-10) )  :min (size (imtestl, 1) ,  (i  +  10) ) 
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max (1, (round (posx (i, j ) ) - 

10) ) : min (size (imtestl, 2) , (round (posx(i,j))+10))),0.9)  ==  1) 
tpxsz(i)  =  tpxsz(i)  +  1; 
tposx (i, tpxsz (i) )  =  posx(i,j); 

end 

end 

end 

%  Then  valleys 
tposx2  =  [ ] ; 

tpx2sz  =  zeros (1, size (imtestl, 1) ) ; 
for  i  =  1 : size (imtestl, 1) 
for  j  =  l:px2sz(i) 

if  (IsLine (imtestl (max (1, (i-10) ) :min (size (imtestl, l),(itl0)),  .. 

max  (1,  (round (posx2  ( i , j ) ) - 

10) ) : min (size (imtestl, 2) , (round (posx2(i,j))+10))),0.9)  ==  1) 
tpx2sz(i)  =  tpx2sz(i)  +  1; 
tposx2 (i, tpx2sz (i) )  =  posx2(i,j); 

end 

end 

end 

disp ( ' Done ' ) ; 
toe 

%  Converting  separate  points  into  connected  line 
disp (' Collecting  into  domain  boundary  lines'); 

[numlines  lines  linesz]  =  FindLines (imtestl, tpx2sz, tposx2) ; 

disp ( ' Done ' ) ; 

toe 

starty  =  1; 

%  Show  filtered  image  with  lines 
figure; 

imshow (umperpx. * (1 : size (iml, 2) ) , umperpx. * (1 : size  (iml, 1) ) , iml ) ; 
xlabel (' Distance  (\mum)'); 
ylabel (' Distance  (\mum)'); 
hold  on; 

for  i  =  1: numlines 

X  =  lines (i, 1 : linesz  (i) , 1)  ; 
y  =  lines (i, 1 : linesz  (i) , 2)  ; 

%  Linear  interpolation 

lines2(i,l: (size (iml, l)-starty+l))  =  in terpl(x,y, start yrsize (imtestl, 1) ) 
plot (umperpx. *lines2 (i, 1 : (size (iml, 1) - 
starty+1) ) , umperpx. * (starty: size (iml, 1) )  ,  ' r :  ' ) ; 
end 
toe 

%  Save  all  information  about  image 

save ( fn ( irow,  i col, 1:12),  'lines',  'lines2',  'linesz',  ' tposx2 ' ,  ' tpx2sz  '  ,  ' tposx ' , 
tpxsz ' ) ; 

end 
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Listing  2:  Fiiter.m 

%  Performs  image  processing  filters  and  derivates 
%  I  =  image 

%  Lpp  and  Lqq  are  calculated  as  eigenvalues  of 
%  Hessian,  [Lxx  Lxy;  Lxy  Lyy] ,  at  each  point 

a 

0 

%  The  image  has  a  Gaussian  filter  applied,  and  then 
%  several  partial  derivatives  (in  x  and  y)  are  calculated. 

%  The  partial  derivatives  are  then  aligned  to  the  principal 
%  curvatures.  Lpp  and  Lqq  are  found  from  the  eigenvalues  of 
%  the  Hessian  matrix.  Finally,  the  blob  function  is  calculated. 

function  [L  Lp  Lq  Lpp  Lqq  blob]  =  Filter (I) 

%  Perform  Gaussian  filter  and  regular  partial  derivatives 
L  =  GaussFilter ( I , 3 , 2 ) ; 

Lx  =  imfilter (double (L)  ,  [-1/12  2/3  0  -2/3  1/12] ,' replicate conv' ) ; 

Ly  =  imfilter (double (L) , [-1/12;  2/3;  0;  -2/3;  1 / 12 ],' replicate conv ') ; 

Lxx  =  imfilter (double (L),[-l/12  4/3  -5/2  4/3  -1/12],' replicate ' , ' conv' ) ; 

Lyy  =  imfilter (double (L) , [-1/12;  4/3;  -5/2;  4/3;  -1/12] ,' replicate ',' conv' ) 
Lxy  =  imfilter (double (L)  ,  [1/4  0  -1/4;  0  0  0;  -1/4  0  1/4]  ,' replicate ',' conv' 
%  Wiener  filter  the  result 
Lx  =  wiener2 (Lx,  [5  5] )  ; 

Ly  =  wiener2 (Ly, [5  5] ) ; 

Lxx  =  wiener2 (Lxx,  [5  5] )  ; 

Lyy  =  wiener2 (Lyy, [5  5] ) ; 

Lxy  =  wiener2 (Lxy,  [5  5] ) ; 

%  Align  to  principal  curvatures 

cosb  =  sqrt (1/2 . * (1+ (Lxx-Lyy)  . / sqrt ( (Lxx-Lyy)  . ^2  +  4 . *Lxy. ^2) ) ) ; 
sinb  =  sqrt  (1/2  .*  (1- (Lxx-Lyy)  . /sqrt  (  (Lxx-Lyy)  .  ^^2  +  4  .  *Lxy.  ^2)  ))  ; 

Lp  =  sinb . *Lx-cosb . *Ly; 

Lq  =  cosb . *Lx+sinb . *Ly; 

%  Second  derivatives  found  from  eigenvalues  of  given  matrix 
Lpp  =  zeros (size (L, 1) , size  (L, 2) )  ; 

Lqq  =  zeros (size (L, 1) , size (L, 2) ) ; 
for  i  =  l:size(L,l) 

for  j  =  1 : size (L, 2 ) 

e  =  eig ( [Lxx (i, j )  Lxy(i,j);  Lxy(i,j)  Lyy(i,j)]); 

Lpp ( i , j )  =  e ( 1 ) ; 

Lqq (i, j  )  =  e (2) ; 

end 

end 

Lpp  =  wiener2 (Lpp,  [5  5] ) ; 

Lqq  =  wiener2 (Lqq, [5  5] ) ; 

%  Calculate  blob  function,  excluding  low  2nd  derivatives 
blob  =  wiener 2 ( (Lpp  >  0.2).* (Lqq  >  0.2).* (Lpp . *Lqq) , [ 5  5 ] ) ; 
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Listing  2.1:  Gaussian. m 

%  Discrete  2-D  Gaussian  function 

%  X  =  X  coordinate 

%  Y  =  Y  coordinate 

%  t  =  t  parameter 

function  y  =  Gaussian (x, y, t) 

y  =  1/ (2  .  *pi  .  *t)  .  *exp  (- (x. ''2+y.  ■^2)  /  (2  .  *t)  )  ; 


Listing  2.2:  GaussFilter.m 

%  Applies  Gaussian  image  filter 
%  I  =  image 

%  m  =  Size  of  filter  on  a  side  divided  by  2 
%  t  =  t  parameter 

g, 

o 

%  Creates  the  specified  Gaussian  filter  and  then 
%  applies  it  to  the  given  image, 
function  M  =  GaussFilter ( I , m, t) 

%  Construct  filter 
for  i  =  -m:m 

for  j  =  -m:m 

f ilterl ( i+m+1 , j +m+l )  =  Gaussian (i, j , t) ; 

end 

end 

%  Use  it 

M  =  imf ilter ( I , filter 1 , ' replicate ' , ' conv ' ) ; 
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Listing  3:  FindBoundary.m 

%  Finds  potential  boundary  points 
%  W  =  window  to  look  for  max  value 
%  Th  =  first  derivative  threshold 
%  Th2  =  second  derivative  threshold 
%  Thb  =  blob  threshold 

g, 

o 

%  For  each  row  in  the  image,  the  potential  ridge  and  valley 
%  points  are  located  using  the  partial  derivatives  and  the 
%  threshold  values.  The  positions  are  then  compiled  into 
%  a  list. 

function  [pos  pos2  psz  p2sz]  = 

FindBoundary ( imdp, imdq, imdpp, imdqq, blob, w, th, th2 , thb) 

%  Search  one  row  at  a  time  in  the  x  direction 
pos  =  zeros (size (imdp, 1) , 200) ; 
pos2  =  zeros (size (imdp, 1) , 200) ; 
for  i  =  1 : size (imdp, 1) 

[p  p2]  = 

FindZeros  ( imdp  ( i ,  :  )  , imdq ( i ,  :  )  , imdpp  ( i ,  :  )  ,  imdqq ( i ,  : ) , blob  ( i ,  :  )  ,  w, th, th2 , thb) 
psz(i)  =  size(p,2); 
p2sz(i)  =  size(p2,2); 
if  (size (p, 1 )  ~=  0 ) 

pos (i, 1 :psz  (i) )  =  p; 

end 

if  (size(p2,l)  ~=  0) 

pos2 (i, 1 :p2sz  (i) )  =  p2; 

end 

end 


Listing  3.1:  Findzeros.m 

%  Finds  zero  crossing  in  the  first  derivatives  and  matches  them 
%  2nd  derivatives  that  meet  the  conditions  for  a  ridge  or 
%  valley 

o 

0 

%  p,q,pp,qq  =  partial  derivatives  of  one  row  of  image 
%  w  =  window  to  look  for  max  value 
%  Th  =  first  derivative  threshold 
%  Th2  =  second  derivative  threshold 
%  Thb  =  blob  threshold 

g, 

o 

%  First,  the  row  of  pixels  is  searched  from  left  to  right  until 
%  a  zero  crossing  in  the  first  derivative  is  found.  Then  the 
%  thresholds  are  compared  to  the  second  derivatives  and  the 
%  blob  function.  If  all  of  them  succeed,  the  criteria  for 
%  ridges  and  valleys  are  applied  and  the  positions  that 
%  satisfy  the  criteria  are  saved. 

function  [pos  pos2]  =  FindZeros (p, q, pp, qq, blob, w, th, th2 , thb) 
pos  =  [] ; 
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pos2  =  [  ] ; 
c  z  =  0 ; 
c  z  2  =  0 ; 
c  =  1; 

while  (c  <  size(p,2)) 

%  Look  for  a  zero  crossing 

while  ( (c  <  size(p,2))  &&  (p(c)*p(c+l)  >  0)  &&  (q(c)*q(c+l)  > 

c  =  c  +  1 ; 

end 

%  When  you  find  one,  match  it  against  the  critera 
if  (c  <  size (p, 2 ) ) 

bmax  =  max (abs (blob (max (1, c-w) :min (size (blob, 2) , c+w) ))) ; 
ppmax  =  max (abs (pp (c : (c+1 ) ) ) ) ; 
qqmax  =  max (abs (qq (c : (c+1 ) ) ) ) ; 

%  Check  if  it's  greater  than  the  threshold  values 
if  ( (bmax  <  thb)  &&  ( (ppmax  >  th2)  | I  (qqmax  >  th2) ) ) 

pzero  =  p (c) *p  (c+1 ) ; 
qzero  =  q (c) *q (c+1 ) ; 
ppmin  =  min (pp (c : (c+1 ) ) ) ; 

[pmaxl  pmaxli]  =  max (abs (p (max ( 1 , c-w) : c) ) ) ; 

[pmaxr  pmaxri]  =  max (abs (p (c+1 :min (size (p, 2) , c+w) ) ) ) ; 
[qmaxl  qmaxli]  =  max (abs (q (max ( 1 , c-w) : c) ) ) ; 

[qmaxr  qmaxri]  =  max (abs (q (c+1 :min (size (p, 2) , c+w) ) ) ) ; 

%  Check  for  ridge 

if  ((pzero  <  0)  &&  (ppmin  <=  0)  &&  (ppmax  >=  qqmax)  && 

th2 )  &&  ... 

(pmaxl  >  th)  &&  (pmaxr  >  th)  &&  (p (max ( 1 , c-w) - 
1+pmaxli) . *p (c+pmaxri)  <  0)) 
cz  =  cz  +  1; 

pos(cz)  =  c-p (c) . / (p (c+1) -p (c) ) ; 

end 

%  Check  for  valley 

if  ((qzero  <  0)  &&  (qqmax  >=  0)  &&  (qqmax  >=  ppmax)  && 

th2 )  &&  ... 

(qmaxl  >  th)  &&  (qmaxr  >  th)  &&  (q (max (1, c-w) - 
1+qmaxli) . *q (c+qmaxri)  <  0)) 
cz2  =  cz2  +  1; 

pos2(cz2)  =  c-q (c) . / (q (c+1 ) -q (c) ) ; 

end 

end 

end 

c  =  c  +  1 ; 

end 


(ppmax  > 


(qqmax  > 
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Listing  4:  IsLine.m 

%  Determines  if  the  pixel  in  the  middle  of  the  image  subsection 
%  is  likely  a  line 

g, 

o 

%  im  =  image  subsection 

%  thp  =  threshold  factor  for  mode  of  image 

g, 

o 

%  First  a  threshold  image  is  created  using  a  fraction  (thp)  of 
%  the  statistical  mode  as  a  cutoff  point.  Starting  in  the  upper 
%  left,  the  threshold  image  is  searched  for  sequences  above 
%  threshold  around  the  border  starting  with  the  top,  then  moving 
%  to  the  right,  bottom,  and  finally  the  left.  If  a  sequence 
%  of  at  least  4  pixels  is  found,  it  is  considered  a  line  exiting 
%  the  subsection.  After  the  search  is  completed,  finishing  at 
%  the  upper  left  again,  the  number  of  sequences  is  counted. 

%  If  there  are  2  sequences,  the  sequence  lengths  are  less  than  11 
%  pixels,  and  the  total  pixels  in  the  subsection  above  threshold 
%  is  less  than  55%  of  the  subsection,  then  the  middle  of  the 
%  subsection  is  likely  to  be  a  part  of  a  line. 

function  A  =  IsLine (im, thp) 

cnt  =  0; 
went  =  0 ; 

%  Create  threshold  image  using  thp*mode(im)  as  divider 
th  =  thp . *mode (mode (double ( im) )) ; 
last  =  double ( im ( 1 , 1 ) )  <=  th; 
wpos  =  [ ] ; 

%  Determines  line  has  started  in  the  upper  left 
if  (last) 

went  =  1 ; 
wid(wcnt)  =  1; 
wpos  =  [ 1 ;  1 ]  ; 

end 

%  Start  from  top  border,  count  the  number  of  the  consecutive 
%  pixels  in  each  line  segment  around  the  border  and  what  size 
%  it  is.  If  it  is  less  than  4  pixels  wide,  disregard 
for  i  =  2:size(im,2) 
if  (last) 

wid(wcnt)  =  wid(wcnt)  +  1; 

end 

if  ( (double ( im ( 1 , i ) )  <=  th)  ~=  last) 
last  =  (double ( im ( 1 , i ) )  <=  th) ; 
if  (last) 

went  =  went  +  1 ; 
wid(wcnt)  =  1 ; 
wpos  =  [1;  i] ; 

else 

if  (wid(wcnt)  <  4) 
went  =  went  -  1; 

else 

mpos  =  ([1;  i ] +wpos ) . /2 ; 
mx  = 

round (1 inspace (round (size(im,l) ./2) ,mpos (2) , size (im, 1) ) ) ; 
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my  = 

round (1 inspace (round (size(im,l) ./2) ,mpos (l),size(im,l))) 
for  k  =  l:size(im,l) 

mim(k)  =  im (my ( k) , mx ( k) )  <  th; 

end 

if  (min (mim)  ~=  i) 

went  =  went  -  1 ; 

end 

end 

end 

ent  =  ent  +  i ; 

end 

end 

%  Do  the  same  on  the  right  border 
for  i  =  i:size(im,i) 
if  (iast) 

wid(wcnt)  =  wid(wcnt)  +  i ; 

end 

if  ( (doubie ( im ( i , size ( im, 2 ) ) )  <=  th)  ~=  iast) 
iast  =  (doubie ( im ( i , size ( im, 2 )) )  <=  th) ; 
if  (iast) 

went  =  went  +  1; 

wid(wcnt)  =  i ; 

wpos  =  [i;  size (im, 2) ] ; 

eise 

if  (wid(wcnt)  <  4) 
went  =  went  -  1 ; 

eise 

mpos  =  ( [i;  size (im, 2) ] +wpos) . /2; 
mx  = 

round (i in space (round (size(im,l) ./2) ,mpos (2),size(im,i))) 
my  = 

round (i inspace (round (size(im,l) ./2) ,mpos (l),size(im,i))) 
for  k  =  i:size(im,i) 

mim(k)  =  im (my ( k) , mx ( k) )  <  th; 

end 

if  (min (mim)  ~=  i) 

went  =  went  -  1 ; 

end 

end 

end 

ent  =  ent  t  1; 

end 

end 

%  And  the  bottom  border 
for  i  =  size(im,2) :-l:i 
if  (iast) 

wid(wcnt)  =  wid(wcnt)  +  i ; 

end 

if  (  (doubie ( im ( size ( im, 1 ), i ) )  <=  th)  ~=  iast) 
iast  =  (doubie ( im ( size ( im, i ), i ) )  <=  th) ; 
if  (iast) 

went  =  went  +  1; 

wid(wcnt)  =  i ; 

wpos  =  [size (im, i) ;  i]; 

eise 

if  (wid(wcnt)  <  4) 
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went 


went 


1; 


else 

mpos  =  ( [size (im, 1) ;  i ] +wpos ) . /2 ; 
mx  = 

round (linspaee (round (size(im,l) ./2) ,mpos (2),size(im,l))) 
my  = 

round (linspaee (round (size(im,l) ./2) ,mpos (l),size(im,l))) 
for  k  =  l:size(im,l) 

mim(k)  =  im (my ( k) , mx ( k) )  <  th; 

end 

if  (min (mim)  ~=  1) 

went  =  went  -  1 ; 

end 

end 

end 

ent  =  ent  +  1; 

end 

end 

%  And  the  left  border 
for  i  =  size(im,l) :-l:l 
if  (last) 

wid(went)  =  wid(went)  +  1; 

end 

if  ( (double ( im ( i , 1 ) )  <=  th)  ~=  last) 
last  =  (double ( im ( i , 1 ) )  <=  th) ; 
if  (last) 

went  =  went  +  1; 
wid(went)  =  1; 
wpos  =  [i;  1]; 

else 

if  (wid(went)  <  4) 
went  =  went  -  1; 

else 

mpos  =  ( [i;  1 ] twpos ) . /2 ; 
mx  = 

round (linspaee (round (size(im,l) .12) ,mpos (2),size(im,l))) 
my  = 

round (linspaee (round (size(im,l) .12)  ,mpos (l),size(im,l))) 
for  k  =  l:size(im,l) 

mim(k)  =  im (my ( k) , mx ( k) )  <  th; 

end 

if  (min (mim)  ~=  1) 

went  =  went  -  1; 

end 

end 

end 

ent  =  ent  t  1; 

end 

end 

%  Cheek  if  it  eonneets  baek  around  to  the  upper  left 
if  ( (double ( im ( 1 , 1 ) )  <=  th)  ~=  last) 
ent  =  ent  +  1; 
elseif  (last) 

wid(l)  =  wid(l)  +  wid(went); 
went  =  went  -  1; 

end 
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o\°  o\°  o\° 


If  there  are  2  line  segments  on  the  border  with  widths  less 
than  11,  and  the  structure  does  not  take  up  more  than  55%  of 
the  image  subsection,  then  it  is  probably  a  line 
A  =  0; 

if  ((went  ==  2)  &&  (max (wid)  <  11)  &&  ( sum ( sum (double ( im)  <= 

th)  )  . / (size (im, 1)  . *size (im, 2) )  <=  0.55)) 

A  =  1; 

end 
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Listing  5:  FindLines.m 

%  Constructs  lines  out  of  separate  boundary  points 

o 

0 

%  imtestl  =  image 

%  tpx2sz  =  size  of  points  making  up  valley 
%  tposx2  =  valley  points 
%  lines  =  constructed  lines 
%  linesz  =  size  of  each  line 
%  numlines  =  number  of  lines 

o 

0 

%  A  cluster  of  every  three  columns  in  the  image  created  and  the 
%  number  of  points  within  each  cluster  is  counted.  Ignoring 
%  isolated  small  clusters,  the  number  of  lines  and  the  location 
%  of  the  clusters  is  found  by  looking  for  groups  of  clusters  with 
%  points  surrounded  by  clusters  with  little  or  no  points.  The 
%  locations  are  then  searched  for  a  point  within  the  cluster, 

%  and  these  points  form  the  beginning  of  the  lines.  Above  and 
%  below  the  points  are  searched  for  other  points  within 
%  8  pixels  horizontally  and  these  new  points  are  connected 
%  to  the  lines.  The  next  points  are  searched  for  within  8 
%  pixels  horizontally  of  the  new  points.  Once  the  top  and  bottom 
%  of  the  image  are  reached,  the  line  is  organized  to  contain  the 
%  points  from  the  top  to  the  bottom  of  the  image. 

function  [numlines  lines  linesz]  =  FindLines (imtestl, tpx2sz, tposx2) 

%  Looks  in  bins  of  three  columns  for  clusters  of  points 
divl  =  3; 

pent  =  zeros (1, ceil (size (imtestl, 2) . /divl)  ) ; 

for  i  =  (round (size (imtestl, 1)  . /2) -100)  :  (round (size (imtestl, 1)  . /2) +100) 
if  (tpx2sz  (i)  >  0) 

for  j  =  l:tpx2sz(i) 

pent (ceil (tposx2 (i, j ). /divl) )  =  pent (ceil (tposx2 (i, j ). /divl) )  + 

1; 

end 

end 

end 

%  Ignore  isolated  points 
pent  =  pent . * (pcnt>7 )  ; 

%  Collect  into  clusters  and  mark  position 
numlines  =  0; 
i  =  1; 

while  (i  <  size (pent, 2 ) ) 

while  ((pcnt(i)  ==  0)  &&  (i  <  size  (pent, 2 )) ) 

i  =  i  +  1  ; 

end 

if  (i  <  size (pent, 2 ) ) 

numlines  =  numlines  +  1; 
linep (numlines )  =  i; 

end 

while  ((pcnt(i)  ~=  0)  &&  (i  <  size  (pent, 2 )) ) 

i  =  i  +  1  ; 

end 
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end 


%  Use  cluster  position  to  find  a  starting  point  for  each  line 
%  looking  somewhere  in  the  center  of  the  image 
lines  =  zeros (numlines, size (imtestl,  1)  ,  2)  ; 
linesz  =  ones ( 1 , numlines ) ; 
for  i  =  1: numlines 
potent  =  0; 
potline  =  [ ] ; 

for  j  =  (round (size (imtestl, 1) . /2) -100) : (round (size (imtestl, 1) . /2) +100) 
for  k  =  l:tpx2sz(j) 

if  (ceil (tposx2 ( j , k) . /divl)  ==  linep(i)) 
potent  =  potent  +  1; 
potline (potent, 1 : 2 )  =  [j;  k] ; 

end 

end 

end 

ment  =  0; 
for  j  =  l:potcnt 
ment  =  0; 
for  k  =  Irpotcnt 

if  (abs ( tposx2 (potline (j,l) ,potline(j,2) )- 
tposx2 (potline ( k, 1 ), potline ( k,  2 )) )  <  4) 
ment  =  ment  +  1; 

end 

if  (ment  >  10) 

potstart  =  j ; 
break; 

end 

end 

if  (ment  >  10) 
break; 

end 

end 

%  Once  you  found  the  point,  go  and  down  for  point  and 
%  connect  points  within  8  pixels  of  each  other 
lineup  =  [potline (potstart, 1 ) ; 
tposx2 (potline (potstart, 1) , potline (potstart,  2) ) ]  ; 
linedown  =  lineup; 

for  j  =  potline (potstart, 1) -1 : -1 : 1 
for  k  =  l:tpx2sz(j) 

if  (tpx2sz ( j )  >0) 

if  (abs (lineup (2, 1) -tposx2 (j , k) )  <  8) 

lineup ( 1 : 2 , 1 :( size ( lineup, 2 ) +1 ) )  =  cat(2,[j; 
tposx2 ( j , k) ] , lineup) ; 

break; 

end 

if  (abs (mean (lineup (2, 1 :min (15, size (lineup, 2) ))) -tposx2 (j , k) )  < 


lineup ( 1 : 2 , 1 :( size ( lineup, 2 ) +1 ) )  =  cat(2,[j; 
tposx2 ( j , k) ] , lineup)  ; 

break; 

end 

end 

if  (tposx2(j,k)  >  lineup (2 , 1 ) +2 ) 
break; 

end 
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end 


end 

for  j  =  (potline (potstart, 1 ) +1 ): size (tposx2 , 1 ) 
for  k  =  l:tpx2sz(j) 

if  (tpx2sz ( j )  >0) 

if  (abs (linedown (2, size (linedown, 2) ) -tposx2 ( j , k) )  <  8) 

linedown (1:2,1: (size (linedown, 2)+l))  =  cat (2, linedown, [ j 

tposx2 ( j , k) ] ) ; 

break; 

end 

if  (abs (mean (linedown (2, max (1, (size(linedown,2)- 
15)  ): size (linedown, 2) )) -tposx2 (j , k) )  <  8) 

linedown (1:2,1: (size (linedown, 2)+l))  =  cat (2, linedown, [ j 

tposx2 ( j , k) ] ) ; 

break; 

end 

end 

if  (tposx2(j,k)  >  linedown (2, size (linedown, 2) ) +2) 
break; 

end 

end 

end 

%  Put  lines  together  (top  and  bottom) 

linesz(i)  =  size ( [lineup (1, : )  linedown (1, 2 : size (linedown, 2) )], 2) +2; 
if  (lineup(l,l)  ~=  1) 

if  (linedown (1, size (linedown, 2) )  ~=  size (imtestl) ) 

lines  (i, 1 : linesz (i) , 1)  =  [1  lineup(l,:) 
linedown (1, 2 : size (linedown,  2) )  size (imtestl, 1) ] ; 

lines (i,  1 : linesz (i)  ,  2)  =  [lineup(2,l)  lineup(2,:) 
linedown (2,2: size (linedown, 2) )  linedown (2, size (linedown, 2) ) ] ; 
else 

linesz (i)  =  linesz(i)-l; 

lines (i, 1 : linesz (i) , 1)  =  [1  lineup(l,:) 
iinedown (1,2: size (linedown,  2) ) ]  ; 

lines (i, 1 : linesz (i) , 2)  =  [lineup(2,l)  lineup(2,:) 
linedown (2,2: size (linedown, 2) ) ] ; 
end 

else 

if  (linedown (1,  size (linedown,  2) )  ~=  size (imtestl) ) 

linesz (i)  =  linesz (i)-l; 
lines  (i, 1 : linesz (i) , 1)  =  [lineup(l,:) 
linedown (1, 2 : size (linedown,  2) )  size (imtestl, 1) ] ; 

lines (i, 1 : linesz  (i) , 2)  =  [lineup(2,:) 
iinedown (2,2: size (linedown, 2) )  linedown (2, size (linedown, 2) ) ] ; 
else 

linesz  (i)  =  linesz  (i) -2; 
lines (i,  1 : linesz (i)  ,  1)  =  [lineup(l,:) 
linedown (1,2: size (linedown,  2) ) ]  ; 

lines (i, 1 : linesz (i) , 2)  =  [lineup(2,:) 
linedown (2,2: size (linedown,  2) ) ]  ; 
end 

end 

end 
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Listing  6:  MaskStitch.m 

%  Collects  the  domain  boundary  lines  from  all  the  images 
%  within  a  superimage  row  and  puts  them  into  a  global 
%  coordinate  system 

o 

0 

%  After  gathering  the  images,  the  registration  information  is  loaded 
%  for  each  image,  and  put  into  a  translation  table,  indicating  how 
%  each  image  needs  to  be  translated  to  get  from  the  first  image  to 
%  the  desire  image's  location.  Then,  after  cleaning  up  the  line  data, 

%  the  line  data  is  translated  to  the  global  coordinate  system  using 
%  the  registration  data.  The  lines  are  then  interpolated  to  line 
%  up  with  integer  vertical  pixels.  Ignore  lines  touching  the  border 
%  of  the  images,  each  new  line  is  added  to  a  master  list  with 
%  the  new  coordinate  system.  As  each  line  is  added,  up  to  the  previous 
%  40  lines  are  compared  to  it,  and  if  the  mean  difference  between 
%  corresponding  points  is  less  than  15,  then  the  lines  are  considered 

%  the  same,  and  the  new  line  is  not  added.  After  all  lines  are  in  the 

%  new  list,  the  range  of  rows  in  the  global  coordinate  system  that 

%  the  lines  are  contained  within  is  determined,  and  then  a  domain 

%  map  is  created.  Using  the  domain  map,  the  FFT  of  each  row  is 
%  calculated  and  the  peak  value  is  saved. 

%  Initialize  parameters 
umperpx  =  0.14244. *2; 
starty  =  1; 
rows  =  3 ; 
cols  =  38; 

%  Get  filenames  of  images 
a  =  dir ( ' * . bmp ' ) ; 
for  cnt  =  l:size(a,l) 

fn (floor ( (cnt-1) /cols) +1, mod (cnt-1, cols) +1, : )  =  a (cnt, 1 ) . name; 

end 

%  Load  registration  data  in  structure,  T,  that  can  be  used 
load  trans 

T  =  zeros (rows, cols, 2)  ; 
for  n  =  1: (rows-1) 

for  m  =  1:  (cols-1) 


T (n, m+1 , 1 )  = 

-trans (n, m, 1 ) 

. tdata . T ( 3 , 2 ) 

1  +T  (n,  m,  1 ) 

-1; 

T (n, m+1 , 2 )  = 

-trans (n, m, 1 ) 

. tdata . T ( 3 , 1 ) 

1  +T  (n,  m,  2  ) 

-1; 

end 

if  (n  ~=  rows-1) 

T(n+l,l,l)  = 

-trans  (n, 1,2) 

. tdata . T ( 3 , 2 ) 

1  +T (n, 1,  1) 

-1; 

T(n+l,l,2)  = 

-trans (n, 1 , 2 ) 

. tdata . T ( 3 , 1 ) 

i+T(n,l,2) 

-1; 

end 

end 

for  m  =  1:  (cols-1) 

T(rows,m,l)  =  - trans (n, m, 2 ) . tdata . T ( 3 , 2 ) +T (rows-1 , m, 1 ) -1 ; 
T(rows,m,2)  =  - trans (n, m, 2 ). tdata . T ( 3 , 1 ) +T (rows-1 , m, 2 ) -1 ; 

end 

T (rows, cols, 1)  =  -Tf . tdata . T ( 3 , 2 ) +T (rows , cols-1 , 1 ) -1 ; 

T (rows, cols, 2)  =  -Tf . tdata . T ( 3 , 1 ) +T (rows , cols-1 , 2 ) -1 ; 

%  For  a  row 
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for  mrow  =  3:3 

%  Mark  top  and  bottom  of  rows 
Trow  =  1+max (T ( 1 ,  : , 1 ) )  ; 

Brow  =  size (iml, 1) +min (T (rows, :,!)); 

T12  =  zeros (1835, Brow) ; 

T12s  =  0; 

dsz  =  zeros (rows, cols) ; 
for  n  =  mrow: mrow; 
for  m  =  1 : cols 

%  Get  information  on  image (n,m) 
load(fn(n,m,l:12) ) ; 
numlines  =  size (lines, 1) ; 

%  Eliminate  extra  points  (multiple 


position 


for  i  =  l:numlines 

if  (sum ( (diff ( lines ( i ,:, 1 ) )  <= 


==  0) ) ) 


points  at  same  vertical 

0)  - (lines (i, 2 : size (lines, 2) ,  1) 


f  =  find (diff ( lines  ( i ,:, 1 ) )  <=  0); 
lines  (i,  (f(l)+l)  :size (lines, 2) , 1)  = 
zeros  (1,  size  (lines,  2)  -  (f  (1)  +1)  +1)  ; 
end 

linesz(i)  =  min ( find ( lines ( i ,:, 1 )  ==  0))-l; 


end 

%  Transfer  local  (image)  line  data  to  global  (superimage) 

%  line  data 

for  i  =  l:numlines 

X  =  lines (i, 1 : linesz  (i) , 1)  ; 
y  =  lines (i, 1 : linesz  (i) , 2)  ; 

Is  =  linesz  (i) ; 

%  Ignore  lines  on  right  or  left  image  border  in 
%  case  they  extend  between  images 
if  ((min(y)  >=  2)  &&  (max(y)  <=  size (iml, 2) -2) ) 

lines2(i,l: (size (iml, l)-starty+l))  = 
interpl (x,y,starty:size (iml,  1)  )  ; 

1  = 


interpl ( (1+T (n,m, 1) ) : (size (iml, 1)+T(n,m, 1) ) ,lines2 (i, :)+T(n,m,2) , . . . 


(ceil(l+T(n,m,l) ) ) :floor(size( iml , 1 ) +T (n, m, 1 ) ) , 'linear ' ) ; 

skip  =  false; 

%  Look  at  up  to  the  last  40  lines  added  and  compare 
%  to  new  line,  if  match  (mean  difference  less  15  pixels 
%  away)  then  ignore  new  line  since  it  was  already 
%  included  from  a  different  image 
for  q  =  max(l, (T12s-40) ) :T12s 

d  =  (l((max(ceil(l+T(n,m,l)),Trange(q,l))- 

ceil ( 1+T (n, m, 1 ) ) +1 ) :  ... 


ceil (1+T (n,m, 1) ) +1) ) 


(min (floor (size (iml, 1)+T(n,m,l)), Trange (q, 2) ) - 


T12 (q,max (ceil (1+T (n,m, 1) ) , Trange (q, 1) ) :min (floor (size (iml, 1)+T(n,m,l)), Trang 
e  (q,  2) )  )  )  ; 

if  (abs (mean (d) )  <=  15) 

dm (n, m, dsz (n, m) +1 )  =  mean(d); 
ds (n, m, dsz (n, m) +1 )  =  std(d); 
dsz(n,m)  =  dsz(n,m)  +  1; 
skip  =  true; 
break; 
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end 


end 

%  If  not  already  added,  include  in  master  list  of  lines 
if  (-skip) 

T12((T12s  +  l),  (ceil(l+T(n,m,l)))  :floor(size (iml, 1)+T(n,m,l)))  =  ... 

1; 

T12s  =  T12s  +  1; 

Trange (T12s, : )  =  [ceil ( 1+T (n, m, 1 ) ) ; 
floor (size (iml, 1)+T(n,m,l))]; 

end 

end 

end 

end 

end 

%  Determine  what  rows  contain  the  lines  and  display  all  the  lines 
figure; 

for  k  =  1 : T12s 
hold  on; 
a  =  1; 
b  =  Brow; 
for  j  =  1 : Brow 

if  (T12(k,j)  ~=  0) 
a  =  j; 

break; 

end 

end 

for  j  =  (a+1) :Brow 

if  (T12(k,j)  ==  0) 
b  =  j-1; 
break; 

end 

end 

plot(T12(k,a:b),a:b, 'b'); 

end 

%  Construct  domain  map  and  record  FFT  peaks  for  each  row 
starty  =  3900; 
endy  =  5250; 

%  starty  =  2250; 

%  endy  =  3500; 

%  starty  =  550; 

%  endy  =  1850; 
looky  =  1300; 
figure; 

for  cy  =  1:50: looky 

mask  =  - 1 .  *ones ( 50 , round (endx+T ( 1 ,  cols ,  2  )  )  )  ; 
for  i  =  1:50 

pos2  =  [startx  T12 ( :  ,  i  +  starty+cy-1 )  '  endx+T ( 1 ,  cols ,  2 )]  ; 
for  j  =  l:2:T12s+l 

mask ( i , round (pos2 ( j ) ) : round (pos2(j+l)))  = 
ones ( 1 , round (pos2 ( j  +1 ) ) -round (pos2 ( j ) ) +1 )  ; 
end 

end 

mask ( 1 : 50 , 1 : startx)  =  zeros ( 50 , startx) ; 

mask(l:50,  (endx+T  ( 1 ,  cols ,  2 ) +1 )  :2.^17)  =  zeros  (50, 2  .  ^'17- 
( endx+T (1, cols, 2) ) )  ; 
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imshow (mask) ; 


image (umperpx . * ( 1 : size (mask, 2 ) ) , umperpx . * ( 1 : size (mask, 1 ) ) , 127 . * (mask+1 ) ) 
xlabel (' Distance  (\mum)'); 
ylabel (' Distance  (\mum)'); 

for  i  =  1:50 

fftl  =  abs (fft (mask (i,  : )  '  ,  2 . ^17)  . /size (mask, 2) ) ; 

[a  b]  =  max  ( f  f  tl  ( 100  :  2  .  17- 100 )  )  ; 
pk (mrow, i+cy-1 )  =  a; 

end 

end 


end 

%  Compare  FFT  peaks  with  ideal  1st  order  QPM 
figure; 

plot(l:size(pk,2) ,pk(3, l:size(pk,2) )  ./  (2. /pi) ) ; 
xlabel (' Verticai  position  (pixel)'); 
ylabel ( ' d  (measured) /d_l ' ) ; 

title (' Actual  d_e_f_f  performance  compared  to  ideal  1st  order  QPM'); 
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Listing  7:  Register.m 

%  Registers  each  image  to  the  images  to  the  right  and  below 
%  then  saves  the  translation  relationships 

g, 

o 

%  startx  =  global  horizontal  position  at  beginning  of  crystal 
%  endx  =  global  horizontal  position  at  end  of  crystal 

g, 

o 

%  After  gathering  all  of  the  image  filenames,  each  image  except 
%  in  the  last  row  or  in  the  last  column  are  loaded  along  with 
%  its  neighbor  to  the  right  and  below.  A  point  on  the  image 
%  and  the  corresponding  point  on  each  of  its  neighbors  is 
%  allowed  to  be  selected.  These  points  are  registered  to  each 
%  other  using  a  cross  correlation  to  improve  precision.  The 
%  image  in  the  last  row,  but  second  last  column  is  then 
%  registered  to  the  image  in  the  last  row  and  column.  Then 
%  a  point  on  the  first  image  is  used  to  show  the  horizontal 
%  position  of  the  beginning  of  the  crystal  and  the  same  with 
%  the  last  column  in  the  first  row  for  the  end  of  the  crystal. 

%  Gathers  filenames  for  each  image 

rows  =  3 ; 

cols  =  38; 

a  =  dir ( ' * . bmp ' ) ; 

for  cnt  =  l:size(a,l) 

fn (floor ( (cnt-1) /cols) +1, mod (cnt-1, cols) +1, : )  =  a (cnt, 1 ) . name; 

end 

%  Load  each  image  and  register  it  to  the  image  to  the  right  and  below 
warning  off; 
colors  =  2; 
overlap  =  0.15; 
iml  =  imread ( f n ( 1 , 1 , : ) ) ; 
iml  =  iml ( : , : , 2 ) ; 
for  rent  =  1: (rows-1) 
for  cent  =  1: (cols-1) 

%  Open  right  and  below  images 
im2  =  imread (fn (rent, ccnt+1, :)) ; 
im3  =  imread (fn (rcnt+1,  cent,  :))  ; 
im2  =  im2 ( : ,  :  ,  2 )  ; 
im3  =  im3 ( : , : , 2 ) ; 

%  Allow  user  to  select  matching  points  in  right  image 
[imlpt, im2pt]  =  cpselect (iml, im2,  'Wait ',  true) ; 

%  Use  correlation  to  improve  registration  precision 
imlptcr  =  epeorr (imlpt, im2pt, iml, im2) ; 

%  Save  data 

ptl2 (rent , cent , 1 , 1 : size ( imlptcr , 1 ), 1 : 2 )  =  imlptcr; 
ptl2 (rent , cent , 2 , 1 : size ( im2pt , 1 ) , 1 : 2 )  =  im2pt; 

%  Convert  data  point  relation  into  transform 
tl2  =  cp2tform (imlptcr, im2pt, ' linear  conformal'); 
trans (rent, cent, 1 )  =  tl2; 

%  Do  again  for  below  image 

[imlpt2, im3pt]  =  cpselect (iml, im3, 'Wait ', true) ; 
imlpt2cr  =  epeorr (imlpt2, im3pt, iml, im3) ; 
ptl3 (rent , cent , 1 , 1 : size ( imlpt2cr, 1 ), 1 : 2 )  =  imlpt2cr; 
ptl3 (rent , cent , 2 , 1 : size ( im3pt , 1 ) , 1 : 2 )  =  im3pt; 
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tl3  =  cp2tform (imlpt2cr, im3pt, ' linear  conformal') 
trans (rent, cent, 2 )  =  tl3; 
trans (rent, cent, 1 ) .tdata.T 
trans (rent, cent, 2 ) .tdata.T 

iml  =  im2 ; 

end 

%  Load  first  image  in  next  superimage  row 
iml  =  imread (fn (rcnt+1, 1,  :  )  )  ; 
iml  =  iml ( : ,  :  ,  2 )  ; 

end 

%  Register  bottom  right  image  with  image  to  left 

iml  =  imread (fn (rows, cols-1,  :))  ; 

iml  =  iml ( : , : , 2 ) ; 

im2  =  imread (fn (rows, cols,  :))  ; 

im2  =  im2 ( : , : , 2 ) ; 

[imlpt, im2pt]  =  cpselect (iml, im2, 'Wait ', true) ; 
imlptcr  =  epeorr (imlpt, im2pt, iml, im2) ; 

Tf  =  cp2tform (imlptcr, im2pt, ' linear  conformal'); 

%  Allow  user  to  select  beginning  of  crystal 
iml  =  imread ( f n ( 1 , 1 , : ) ) ; 
iml  =  iml ( : , : , 2 ) ; 
figure;  imshow(iml); 

[x  y]  =  ginput (1) ; 
startx  =  x; 

%  and  end  of  crystal 

iml  =  imread ( fn ( 1 ,  cols ,:))  ; 

iml  =  iml ( : , : , 2 ) ; 

figure;  imshow(iml); 

[x  y]  =  ginput (1) ; 
endx  =  x; 

save  Trans  ptl2  ptl3  trans  Tf  startx  endx 
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Listing  8:  Review.m 

%  Checks  to  make  sure  that  there  are  no  duplicate  domain  boundary 
%  lines  within  the  master  list.  If  there  is,  they  are  removed 

g, 

o 

%  iml  =  image 

%  numlines  =  number  of  lines  in  image 
%  lines  =  list  of  lines  (collection  of  data  points) 

%  lines2  =  list  of  interpolated  lines  (horizontal  position 
%  for  each  vertical  pixel  in  line) 

o 

0 

%  For  each  line,  the  next  line  is  compared  point  by  point.  Each  point 
%  within  1  pixel  of  the  next  line's  corresponding  point  is  counted. 

%  If  90%  of  the  points  match  by  this  criteria,  the  lines  are 
%  considered  the  same  and  one  of  them  is  eliminated. 

function  [numlines  lines  linesz  lines2]  = 

Review (iml, numlines, lines, linesz, lines2) ; 

%  Search  through  each  line 
for  i  =  1: (numlines-1) 
if  (i  >  numlines-1) 
break; 

end 

matches  =  0; 

%  Check  point  in  line  for  a  match  in  another  line 
for  j  =  1 : size (iml, 1) 

if  (abs (lines2 (i, j ) -lines2 (i+1, j ) )  <=  1) 
matches  =  matches  +  1; 

end 

end 

%  Try  to  fix  mistake  (eliminate  duplicates) 
if  (matches  >  0 . 9*size (iml, 1) ) 
if  (i  <  numlines-1) 

linesz  =  [linesz  (l:i)  linesz (i+2 rnumlines) ] ; 
lines  = 

cat(l,lines(l:i,l:size(lines,2)  ,1:2)  ,lines( (i+2)  : numlines , l:size(lines,2) ,1:2 

) )  ; 

lines2  = 

cat(l,lines2 (l:i, l:size(lines,2) ) ,lines2 ( (i+2) : numlines, l:size(lines,2) ) ) ; 
numlines  =  numlines  -  1; 

else 

linesz  =  [linesz (l:i)  linesz (i+2 :numlines) ] ; 
lines  =  lines ( 1 : i , 1 : size ( lines , 2 ), 1 : 2 ) ; 
lines2  =  lines2 (1 : i, 1 : size (lines, 2) )  ; 
numlines  =  numlines  -  1; 

end 

end 

end 
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